home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / LFIT.DEM < prev    next >
Text File  |  1991-04-29  |  3KB  |  122 lines
  1. PROGRAM d14r2(input,output);
  2. (* driver for routine LFIT *)
  3. CONST
  4.    npt=100;
  5.    spread=0.1;
  6.    nterm=3;   
  7. TYPE
  8.    glcovar = ARRAY [1..nterm,1..nterm] OF real;
  9.    glnpbynp = glcovar;
  10.    glnpbymp = ARRAY [1..nterm,1..1] OF real;
  11.    gllista = ARRAY [1..nterm] OF integer;
  12.    glndata = ARRAY [1..npt] OF real;
  13.    glmma = ARRAY [1..nterm] OF real;
  14. VAR
  15.    gliset : integer;
  16.    glgset : real;
  17.    glinext,glinextp : integer;
  18.    glma : ARRAY [1..55] OF real;
  19.    chisq : real;
  20.    i,ii,idum,j,mfit : integer;
  21.    lista : gllista;
  22.    a : glmma;
  23.    covar : glcovar;
  24.    x,y,sig : glndata;
  25.  
  26. PROCEDURE funcs(x: real; VAR afunc: glmma; mma: integer);
  27. (* Programs using FUNCS must define the type
  28. TYPE
  29.    glmma = ARRAY [1..mma] OF real;
  30. in the main routine. *)
  31. VAR
  32.    i : integer;
  33. BEGIN
  34.    afunc[1] := 1.0;
  35.    FOR i := 2 to mma DO BEGIN
  36.       afunc[i] := x*afunc[i-1]
  37.    END
  38. END;
  39.  
  40. (*$I MODFILE.PAS *)
  41. (*$I RAN3.PAS *)
  42.  
  43. (*$I GASDEV.PAS *)
  44.  
  45. (*$I GAUSSJ.PAS *)
  46.  
  47. (*$I COVSRT.PAS *)
  48.  
  49. (*$I LFIT.PAS *)
  50.  
  51. BEGIN
  52.    gliset := 0;
  53.    idum := -911;
  54.    FOR i := 1 to npt DO BEGIN
  55.       x[i] := 0.1*i;
  56.       y[i] := nterm;
  57.       FOR j := nterm-1 DOWNTO 1 DO BEGIN
  58.          y[i] := j+y[i]*x[i]
  59.       END;
  60.       y[i] := y[i]+spread*gasdev(idum);
  61.       sig[i] := spread
  62.    END;
  63.    mfit := nterm;
  64.    FOR i := 1 to mfit DO BEGIN
  65.       lista[i] := i
  66.    END;
  67.    lfit(x,y,sig,npt,a,nterm,lista,mfit,covar,nterm,chisq);
  68.    writeln;
  69.    writeln('parameter':9,'uncertainty':23);
  70.    FOR i := 1 to nterm DO BEGIN
  71.       writeln('a[':4,i:1,'] = ',a[i]:8:6,sqrt(covar[i,i]):12:6)
  72.    END;
  73.    writeln('chi-squared = ',chisq:12);
  74.    writeln('full covariance matrix');
  75.    FOR i := 1 to nterm DO BEGIN
  76.       FOR j := 1 to nterm DO write(covar[i,j]:12);
  77.       writeln
  78.    END;
  79.    writeln;
  80.    writeln('press RETURN to continue...');
  81.    readln;
  82. (* now test the LISTA feature *)
  83.    FOR i := 1 to nterm DO BEGIN
  84.       lista[i] := nterm+1-i
  85.    END;
  86.    lfit(x,y,sig,npt,a,nterm,lista,mfit,covar,nterm,chisq);
  87.    writeln('parameter':9,'uncertainty':23);
  88.    FOR i := 1 to nterm DO BEGIN
  89.       writeln('a[':4,i:1,'] = ',a[i]:8:6,sqrt(covar[i,i]):12:6)
  90.    END;
  91.    writeln('chi-squared = ',chisq:12);
  92.    writeln('full covariance matrix');
  93.    FOR i := 1 to nterm DO BEGIN
  94.       FOR j := 1 to nterm DO write(covar[i,j]:12);
  95.       writeln
  96.    END;
  97.    writeln;
  98.    writeln('press RETURN to continue...');
  99.    readln;
  100. (* now check results of restricting fit parameters *)
  101.    ii := 1;
  102.    FOR i := 1 to nterm DO BEGIN
  103.       IF ((i MOD 2) = 1) THEN BEGIN
  104.          lista[ii] := i;
  105.          ii := ii+1
  106.       END
  107.    END;
  108.    mfit := ii-1;
  109.    lfit(x,y,sig,npt,a,nterm,lista,mfit,covar,nterm,chisq);
  110.    writeln('parameter':9,'uncertainty':23);
  111.    FOR i := 1 to nterm DO BEGIN
  112.       writeln('a[':4,i:1,'] = ',a[i]:8:6,sqrt(covar[i,i]):12:6)
  113.    END;
  114.    writeln('chi-squared = ',chisq:12);
  115.    writeln('full covariance matrix');
  116.    FOR i := 1 to nterm DO BEGIN
  117.       FOR j := 1 to nterm DO write(covar[i,j]:12);
  118.       writeln
  119.    END;
  120.    writeln
  121. END.
  122.